# Load packages
  library(tidyverse)
  library(broom)
  library(lmtest)
  library(sandwich)
  library(stargazer)
  
# Calculate robust confidence intervals
  se_robust <- function(x)
    coeftest(x, vcov = vcovHC(x, type = "HC0"))[, "Std. Error"]
  p_robust <- function(x)
    coeftest(x, vcov = vcovHC(x, type = "HC0"))[, "Pr(>|t|)"]

  #### MEDW 2014 Belgian National Election Study
  
  # Load data
  load("ObsStudy_00_data_Belgium_2014.RData")
  dt <- x

  # Party rating
  colnames(dt)[seq(51, 58, 1)] <- c("rat_PS", "rat_MR", "rat_CDH", "rat_Ecolo", "rat_PTB", "rat_FDF", "rat_PP", "rat_FN")
  colnames(dt)[seq(291, 298, 1)] <- c("rat_N-VA", "rat_CD&V", "rat_OpenVLD", "rat_SP.A", "rat_VlaamsBelang", "rat_Groen", "rat_LDD", "rat_PVDA")
  
  # Set value 99 (Ne sait pas/Don´t know) missing
  dt[, 51:58][dt[, 51:58] == 99] <- NA
  dt[, 291:298][dt[, 291:298] == 99] <- NA
  
  # Coalition evaluations/preferred coalition
  colnames(dt)[seq(62, 65, 1)] <- c("rat_coal_PS_SP.A_CD&V_CDH_MR_OpenVLD", "rat_coal_N-VA_CD&V_CDH_MR_OpenVLD", "rat_coal_PS_SP.A_CD&V_CDH_MR_OpenVLD_Ecolo_Groen", "rat_coal_CD&V_CDH_PS_SP.A_Ecolo_Groen")
  
  # Set value 99 (Ne sait pas/Don´t know) missing
  dt[, 62:65][dt[, 62:65] == 99] <- NA

  # Coalition likelihood
  colnames(dt)[seq(66, 69, 1)] <- c("coallik_PS_SP.A_CD&V_CDH_MR_OpenVLD", "coallik_N-VA_CD&V_CDH_MR_OpenVLD", "coallik_PS_SP.A_CD&V_CDH_MR_OpenVLD_Ecolo_Groen", "coallik_CD&V_CDH_PS_SP.A_Ecolo_Groen")
  
  # Set value 99 (Ne sait pas/Don´t know) missing
  dt[, 66:69][dt[, 66:69] == 99] <- NA
  
  # Dependent Variable: nominal vote choice
  # Set missing codes  to NA 
  dt[, 22][dt[, 22] >= 88] <- NA
  dt$vote <- dt$Q8A
  
  dt$vote[dt$vote==1] <- "PS"
  dt$vote[dt$vote==2] <- "MR"
  dt$vote[dt$vote==3] <- "CDH"
  dt$vote[dt$vote==4] <- "Ecolo"
  dt$vote[dt$vote==5] <- "PTB"
  dt$vote[dt$vote==6] <- "FDF"
  dt$vote[dt$vote==7] <- "PP"
  dt$vote[dt$vote==8] <- "FN"
  dt$vote[dt$vote==9] <- "N-VA"
  dt$vote[dt$vote==10] <- "CD&V"
  dt$vote[dt$vote==11] <- "OpenVLD"
  dt$vote[dt$vote==12] <- "SP.A"
  dt$vote[dt$vote==13] <- "VlaamsBelang"
  dt$vote[dt$vote==14] <- "Groen"
  dt$vote[dt$vote==15] <- "LDD"
  dt$vote[dt$vote==16] <- "PVDA"
  
  # Generate vote choice columns based on vote variable
  dt$ptv_PS <- ifelse(dt$vote == "PS", 1, 0)
  dt$ptv_SP.A <- ifelse(dt$vote == "SP.A", 1, 0)
  dt$ptv_CDV <- ifelse(dt$vote == "CD&V", 1, 0)
  dt$ptv_CDH <- ifelse(dt$vote == "CDH", 1, 0)
  dt$ptv_MR <- ifelse(dt$vote == "MR", 1, 0)
  dt$ptv_OpenVLD <- ifelse(dt$vote == "OpenVLD", 1, 0)
  dt$ptv_NVA <- ifelse(dt$vote == "N-VA", 1, 0)
  dt$ptv_Ecolo <- ifelse(dt$vote == "Ecolo", 1, 0)
  dt$ptv_Groen <- ifelse(dt$vote == "Groen", 1, 0)
  
  # Coalition evaluation
  Z <- dt %>% select(contains("rat_coal")) %>%
    mutate_all(as.numeric) %>%
    mutate_all(funs(scales::rescale(.,to = c(0, 1))))
  
  rat_coal_PS <- Z %>% select(`rat_coal_PS_SP.A_CD&V_CDH_MR_OpenVLD`, `rat_coal_PS_SP.A_CD&V_CDH_MR_OpenVLD_Ecolo_Groen`, `rat_coal_CD&V_CDH_PS_SP.A_Ecolo_Groen`) %>% as.matrix()
  rat_coal_SP.A <- Z %>% select(`rat_coal_PS_SP.A_CD&V_CDH_MR_OpenVLD`, `rat_coal_PS_SP.A_CD&V_CDH_MR_OpenVLD_Ecolo_Groen`, `rat_coal_CD&V_CDH_PS_SP.A_Ecolo_Groen`) %>% as.matrix()
  rat_coal_CDV <- Z %>% select(`rat_coal_PS_SP.A_CD&V_CDH_MR_OpenVLD`, `rat_coal_N-VA_CD&V_CDH_MR_OpenVLD`, `rat_coal_PS_SP.A_CD&V_CDH_MR_OpenVLD_Ecolo_Groen`, `rat_coal_CD&V_CDH_PS_SP.A_Ecolo_Groen`) %>% as.matrix()
  rat_coal_CDH <- Z %>% select(`rat_coal_PS_SP.A_CD&V_CDH_MR_OpenVLD`, `rat_coal_N-VA_CD&V_CDH_MR_OpenVLD`, `rat_coal_PS_SP.A_CD&V_CDH_MR_OpenVLD_Ecolo_Groen`, `rat_coal_CD&V_CDH_PS_SP.A_Ecolo_Groen`) %>% as.matrix()
  rat_coal_MR <- Z %>% select(`rat_coal_PS_SP.A_CD&V_CDH_MR_OpenVLD`, `rat_coal_N-VA_CD&V_CDH_MR_OpenVLD`, `rat_coal_PS_SP.A_CD&V_CDH_MR_OpenVLD_Ecolo_Groen`) %>% as.matrix()
  rat_coal_OpenVLD <- Z %>% select(`rat_coal_PS_SP.A_CD&V_CDH_MR_OpenVLD`, `rat_coal_N-VA_CD&V_CDH_MR_OpenVLD`, `rat_coal_PS_SP.A_CD&V_CDH_MR_OpenVLD_Ecolo_Groen`) %>% as.matrix()
  rat_coal_Ecolo <- Z %>% select(`rat_coal_PS_SP.A_CD&V_CDH_MR_OpenVLD_Ecolo_Groen`, `rat_coal_CD&V_CDH_PS_SP.A_Ecolo_Groen`) %>% as.matrix()
  rat_coal_Groen <- Z %>% select(`rat_coal_PS_SP.A_CD&V_CDH_MR_OpenVLD_Ecolo_Groen`, `rat_coal_CD&V_CDH_PS_SP.A_Ecolo_Groen`) %>% as.matrix()
  
  
  # Coalition likelihood
  gamma <- dt %>% select(contains("coallik_")) 
  
  
  coal_lik_PS <- gamma %>% select(`coallik_PS_SP.A_CD&V_CDH_MR_OpenVLD`, `coallik_PS_SP.A_CD&V_CDH_MR_OpenVLD_Ecolo_Groen`, `coallik_CD&V_CDH_PS_SP.A_Ecolo_Groen`) %>%
    mutate(sum_exp = exp(`coallik_PS_SP.A_CD&V_CDH_MR_OpenVLD`) + exp(`coallik_PS_SP.A_CD&V_CDH_MR_OpenVLD_Ecolo_Groen`) + exp(`coallik_CD&V_CDH_PS_SP.A_Ecolo_Groen`),
           `coallik_PS_SP.A_CD&V_CDH_MR_OpenVLD` = exp(`coallik_PS_SP.A_CD&V_CDH_MR_OpenVLD`)/sum_exp,
           `coallik_PS_SP.A_CD&V_CDH_MR_OpenVLD_Ecolo_Groen` = exp(`coallik_PS_SP.A_CD&V_CDH_MR_OpenVLD_Ecolo_Groen`)/sum_exp,
           `coallik_CD&V_CDH_PS_SP.A_Ecolo_Groen` = exp(`coallik_CD&V_CDH_PS_SP.A_Ecolo_Groen`)/sum_exp) %>%
   select(-sum_exp) %>%
    as.matrix()
  
  coal_lik_SP.A <- gamma %>% select(`coallik_PS_SP.A_CD&V_CDH_MR_OpenVLD`, `coallik_PS_SP.A_CD&V_CDH_MR_OpenVLD_Ecolo_Groen`, `coallik_CD&V_CDH_PS_SP.A_Ecolo_Groen`) %>%
    mutate(sum_exp = exp(`coallik_PS_SP.A_CD&V_CDH_MR_OpenVLD`) + exp(`coallik_PS_SP.A_CD&V_CDH_MR_OpenVLD_Ecolo_Groen`) + exp(`coallik_CD&V_CDH_PS_SP.A_Ecolo_Groen`),
           `coallik_PS_SP.A_CD&V_CDH_MR_OpenVLD` = exp(`coallik_PS_SP.A_CD&V_CDH_MR_OpenVLD`)/sum_exp,
           `coallik_PS_SP.A_CD&V_CDH_MR_OpenVLD_Ecolo_Groen` = exp(`coallik_PS_SP.A_CD&V_CDH_MR_OpenVLD_Ecolo_Groen`)/sum_exp,
           `coallik_CD&V_CDH_PS_SP.A_Ecolo_Groen` = exp(`coallik_CD&V_CDH_PS_SP.A_Ecolo_Groen`)/sum_exp) %>%
    select(-sum_exp) %>%
    as.matrix()
  
  coal_lik_CDV <- gamma %>% select(`coallik_PS_SP.A_CD&V_CDH_MR_OpenVLD`, `coallik_N-VA_CD&V_CDH_MR_OpenVLD`, `coallik_PS_SP.A_CD&V_CDH_MR_OpenVLD_Ecolo_Groen`, `coallik_CD&V_CDH_PS_SP.A_Ecolo_Groen`) %>%
    mutate(sum_exp = exp(`coallik_PS_SP.A_CD&V_CDH_MR_OpenVLD`) +  exp(`coallik_N-VA_CD&V_CDH_MR_OpenVLD`) + exp(`coallik_PS_SP.A_CD&V_CDH_MR_OpenVLD_Ecolo_Groen`) + exp(`coallik_CD&V_CDH_PS_SP.A_Ecolo_Groen`),
           `coallik_PS_SP.A_CD&V_CDH_MR_OpenVLD` = exp(`coallik_PS_SP.A_CD&V_CDH_MR_OpenVLD`)/sum_exp,
           `coallik_N-VA_CD&V_CDH_MR_OpenVLD` = exp(`coallik_N-VA_CD&V_CDH_MR_OpenVLD`)/sum_exp,
           `coallik_PS_SP.A_CD&V_CDH_MR_OpenVLD_Ecolo_Groen` = exp(`coallik_PS_SP.A_CD&V_CDH_MR_OpenVLD_Ecolo_Groen`)/sum_exp,
           `coallik_CD&V_CDH_PS_SP.A_Ecolo_Groen` = exp(`coallik_CD&V_CDH_PS_SP.A_Ecolo_Groen`)/sum_exp) %>%
    select(-sum_exp) %>%
    as.matrix()
  
  coal_lik_CDH <- gamma %>% select(`coallik_PS_SP.A_CD&V_CDH_MR_OpenVLD`, `coallik_N-VA_CD&V_CDH_MR_OpenVLD`, `coallik_PS_SP.A_CD&V_CDH_MR_OpenVLD_Ecolo_Groen`, `coallik_CD&V_CDH_PS_SP.A_Ecolo_Groen`) %>%
    mutate(sum_exp = exp(`coallik_PS_SP.A_CD&V_CDH_MR_OpenVLD`) +  exp(`coallik_N-VA_CD&V_CDH_MR_OpenVLD`) + exp(`coallik_PS_SP.A_CD&V_CDH_MR_OpenVLD_Ecolo_Groen`) + exp(`coallik_CD&V_CDH_PS_SP.A_Ecolo_Groen`),
           `coallik_PS_SP.A_CD&V_CDH_MR_OpenVLD` = exp(`coallik_PS_SP.A_CD&V_CDH_MR_OpenVLD`)/sum_exp,
           `coallik_N-VA_CD&V_CDH_MR_OpenVLD` = exp(`coallik_N-VA_CD&V_CDH_MR_OpenVLD`)/sum_exp,
           `coallik_PS_SP.A_CD&V_CDH_MR_OpenVLD_Ecolo_Groen` = exp(`coallik_PS_SP.A_CD&V_CDH_MR_OpenVLD_Ecolo_Groen`)/sum_exp,
           `coallik_CD&V_CDH_PS_SP.A_Ecolo_Groen` = exp(`coallik_CD&V_CDH_PS_SP.A_Ecolo_Groen`)/sum_exp) %>%
    select(-sum_exp) %>%
    as.matrix()
  
  coal_lik_MR <- gamma %>% select(`coallik_PS_SP.A_CD&V_CDH_MR_OpenVLD`, `coallik_N-VA_CD&V_CDH_MR_OpenVLD`, `coallik_PS_SP.A_CD&V_CDH_MR_OpenVLD_Ecolo_Groen`) %>%
    mutate(sum_exp = exp(`coallik_PS_SP.A_CD&V_CDH_MR_OpenVLD`) +  exp(`coallik_N-VA_CD&V_CDH_MR_OpenVLD`) + exp(`coallik_PS_SP.A_CD&V_CDH_MR_OpenVLD_Ecolo_Groen`),
           `coallik_PS_SP.A_CD&V_CDH_MR_OpenVLD` = exp(`coallik_PS_SP.A_CD&V_CDH_MR_OpenVLD`)/sum_exp,
           `coallik_N-VA_CD&V_CDH_MR_OpenVLD` = exp(`coallik_N-VA_CD&V_CDH_MR_OpenVLD`)/sum_exp,
           `coallik_PS_SP.A_CD&V_CDH_MR_OpenVLD_Ecolo_Groen` = exp(`coallik_PS_SP.A_CD&V_CDH_MR_OpenVLD_Ecolo_Groen`)/sum_exp) %>%
    select(-sum_exp) %>%
    as.matrix()
  
  coal_lik_OpenVLD <- gamma %>% select(`coallik_PS_SP.A_CD&V_CDH_MR_OpenVLD`, `coallik_N-VA_CD&V_CDH_MR_OpenVLD`, `coallik_PS_SP.A_CD&V_CDH_MR_OpenVLD_Ecolo_Groen`) %>%
    mutate(sum_exp = exp(`coallik_PS_SP.A_CD&V_CDH_MR_OpenVLD`) +  exp(`coallik_N-VA_CD&V_CDH_MR_OpenVLD`) + exp(`coallik_PS_SP.A_CD&V_CDH_MR_OpenVLD_Ecolo_Groen`),
           `coallik_PS_SP.A_CD&V_CDH_MR_OpenVLD` = exp(`coallik_PS_SP.A_CD&V_CDH_MR_OpenVLD`)/sum_exp,
           `coallik_N-VA_CD&V_CDH_MR_OpenVLD` = exp(`coallik_N-VA_CD&V_CDH_MR_OpenVLD`)/sum_exp,
           `coallik_PS_SP.A_CD&V_CDH_MR_OpenVLD_Ecolo_Groen` = exp(`coallik_PS_SP.A_CD&V_CDH_MR_OpenVLD_Ecolo_Groen`)/sum_exp) %>%
    select(-sum_exp) %>%
    as.matrix()
  
  
  coal_lik_Ecolo <- gamma %>% select(`coallik_PS_SP.A_CD&V_CDH_MR_OpenVLD_Ecolo_Groen`, `coallik_CD&V_CDH_PS_SP.A_Ecolo_Groen`) %>%
    mutate(sum_exp = exp(`coallik_PS_SP.A_CD&V_CDH_MR_OpenVLD_Ecolo_Groen`) + exp(`coallik_CD&V_CDH_PS_SP.A_Ecolo_Groen`),
           `coallik_PS_SP.A_CD&V_CDH_MR_OpenVLD_Ecolo_Groen` = exp(`coallik_PS_SP.A_CD&V_CDH_MR_OpenVLD_Ecolo_Groen`)/sum_exp,
           `coallik_CD&V_CDH_PS_SP.A_Ecolo_Groen` = exp(`coallik_CD&V_CDH_PS_SP.A_Ecolo_Groen`)/sum_exp) %>%
    select(-sum_exp) %>%
    as.matrix()
  
  coal_lik_Groen <- gamma %>% select(`coallik_PS_SP.A_CD&V_CDH_MR_OpenVLD_Ecolo_Groen`, `coallik_CD&V_CDH_PS_SP.A_Ecolo_Groen`) %>%
    mutate(sum_exp = exp(`coallik_PS_SP.A_CD&V_CDH_MR_OpenVLD_Ecolo_Groen`) + exp(`coallik_CD&V_CDH_PS_SP.A_Ecolo_Groen`),
           `coallik_PS_SP.A_CD&V_CDH_MR_OpenVLD_Ecolo_Groen` = exp(`coallik_PS_SP.A_CD&V_CDH_MR_OpenVLD_Ecolo_Groen`)/sum_exp,
           `coallik_CD&V_CDH_PS_SP.A_Ecolo_Groen` = exp(`coallik_CD&V_CDH_PS_SP.A_Ecolo_Groen`)/sum_exp) %>%
    select(-sum_exp) %>%
    as.matrix()
  
  
  # Mean-variance model
  EV <- function(gamma,Z){
    gamma %*% Z
  }
  
  Vari <- function(gamma,Z){
    gamma %*% ((Z - as.numeric(EV(gamma,Z)))^2)
  }
  
  N <- nrow(dt)
  
  
  M_PS <- sapply(1:N, function(i) EV(coal_lik_PS[i,], rat_coal_PS[i,]))
  V_PS <- sapply(1:N, function(i) Vari(coal_lik_PS[i,], rat_coal_PS[i,]))
  
  M_SP.A <- sapply(1:N, function(i) EV(coal_lik_SP.A[i,], rat_coal_SP.A[i,]))
  V_SP.A <- sapply(1:N, function(i) Vari(coal_lik_SP.A[i,], rat_coal_SP.A[i,]))
  
  M_CDV <- sapply(1:N, function(i) EV(coal_lik_CDV[i,], rat_coal_CDV[i,]))
  V_CDV <- sapply(1:N, function(i) Vari(coal_lik_CDV[i,], rat_coal_CDV[i,]))
  
  M_CDH <- sapply(1:N, function(i) EV(coal_lik_CDH[i,], rat_coal_CDH[i,]))
  V_CDH <- sapply(1:N, function(i) Vari(coal_lik_CDH[i,], rat_coal_CDH[i,])) 
  
  M_MR <- sapply(1:N, function(i) EV(coal_lik_MR[i,], rat_coal_MR[i,]))
  V_MR <- sapply(1:N, function(i) Vari(coal_lik_MR[i,], rat_coal_MR[i,]))
  
  M_OpenVLD <- sapply(1:N, function(i) EV(coal_lik_OpenVLD[i,], rat_coal_OpenVLD[i,]))
  V_OpenVLD <- sapply(1:N, function(i) Vari(coal_lik_OpenVLD[i,], rat_coal_OpenVLD[i,]))
  
  M_Ecolo <- sapply(1:N, function(i) EV(coal_lik_Ecolo[i,], rat_coal_Ecolo[i,]))
  V_Ecolo <- sapply(1:N, function(i) Vari(coal_lik_Ecolo[i,], rat_coal_Ecolo[i,]))
  
  M_Groen <- sapply(1:N, function(i) EV(coal_lik_Groen[i,], rat_coal_Groen[i,]))
  V_Groen <- sapply(1:N, function(i) Vari(coal_lik_Groen[i,], rat_coal_Groen[i,]))
  
  # Gender
  dt$male <- ifelse(dt$GEND == 1, 1, 0)
  
  dt$Q46A
  
  # PID
  if (!exists("robustnesscheck_pid")) {
    robustnesscheck_pid <- FALSE # PID as robustness?
  }
  dt$pid <- dt$Q46A
  dt$pid[dt$Q46==2] <- 0
  
  dt$pid_PS <- ifelse(dt$pid==1, 1, 0)
  dt$pid_SP.A <- ifelse(dt$pid==12, 1, 0)
  dt$pid_CDV <- ifelse(dt$pid==10, 1, 0)
  dt$pid_CDH <- ifelse(dt$pid==3, 1, 0)
  dt$pid_MR <- ifelse(dt$pid==2, 1, 0)
  dt$pid_OpenVLD <- ifelse(dt$pid==11, 1, 0)
  dt$pid_Ecolo <- ifelse(dt$pid==4, 1, 0)
  dt$pid_Groen <- ifelse(dt$pid==14, 1, 0)

  d <- data.frame(
    "ptv_PS" = dt$ptv_PS,
    "ptv_SP.A" = dt$ptv_SP.A,
    "ptv_CDV" = dt$ptv_CDV,
    "ptv_CDH" = dt$ptv_CDH,
    "ptv_MR" = dt$ptv_MR,
    "ptv_OpenVLD" = dt$ptv_OpenVLD,
    "ptv_Ecolo" = dt$ptv_Ecolo,
    "ptv_Groen" = dt$ptv_Groen,
    "rat.PS" = dt$rat_PS,
    "rat.SP.A" = dt$rat_SP.A,
    "rat.CDV" = dt$`rat_CD&V`,
    "rat.CDH" = dt$rat_CDH,
    "rat.MR" = dt$rat_MR,
    "rat.OpenVLD" = dt$rat_OpenVLD,
    "rat.Ecolo" = dt$rat_Ecolo,
    "rat.Groen" = dt$rat_Groen,
    "pid_PS" = dt$pid_PS,
    "pid_SP.A" = dt$pid_SP.A,
    "pid_CDV" = dt$pid_CDV,
    "pid_CDH" = dt$pid_CDH,
    "pid_MR" = dt$pid_MR,
    "pid_OpenVLD" = dt$pid_OpenVLD,
    "pid_Ecolo" = dt$pid_Ecolo,
    "pid_Groen" = dt$pid_Groen,
    "lotterymean.PS" = M_PS,
    "lotteryvariance.PS" = V_PS,
    "lotterymean.SP.A" = M_SP.A,
    "lotteryvariance.SP.A" = V_SP.A,
    "lotterymean.CDV" = M_CDV,
    "lotteryvariance.CDV" = V_CDV,
    "lotterymean.CDH" = M_CDH,
    "lotteryvariance.CDH" = V_CDH,
    "lotterymean.MR" = M_MR,
    "lotteryvariance.MR" = V_MR,
    "lotterymean.OpenVLD" = M_OpenVLD,
    "lotteryvariance.OpenVLD" = V_OpenVLD,
    "lotterymean.Ecolo" = M_Ecolo,
    "lotteryvariance.Ecolo" = V_Ecolo,
    "lotterymean.Groen" = M_Groen,
    "lotteryvariance.Groen" = V_Groen,
    "sex" = dt$male,
    "age" = dt$AGE,
    "edu" = dt$SD4
  )
  
  d$lotteryvariance.PS <- scales::rescale(d$lotteryvariance.PS, to = c(0, 1), from = c(0,0.25))
  d$lotteryvariance.SP.A <- scales::rescale(d$lotteryvariance.SP.A, to = c(0, 1), from = c(0,0.25))
  d$lotteryvariance.CDV <- scales::rescale(d$lotteryvariance.CDV, to = c(0, 1), from = c(0,0.25))
  d$lotteryvariance.CDH <- scales::rescale(d$lotteryvariance.CDH, to = c(0, 1), from = c(0,0.25))
  d$lotteryvariance.MR <- scales::rescale(d$lotteryvariance.MR, to = c(0, 1), from = c(0,0.25))
  d$lotteryvariance.OpenVLD <- scales::rescale(d$lotteryvariance.OpenVLD, to = c(0, 1), from = c(0,0.25))
  d$lotteryvariance.Ecolo <- scales::rescale(d$lotteryvariance.Ecolo, to = c(0, 1), from = c(0,0.25))
  d$lotteryvariance.Groen <- scales::rescale(d$lotteryvariance.Groen, to = c(0, 1), from = c(0,0.25))

# Save model results
  summary(m1 <- lm(as.formula(paste("ptv_PS ~ lotteryvariance.PS + lotterymean.PS + rat.PS + sex + as.factor(edu) + age", 
                                    if (robustnesscheck_pid) "+ pid_PS" else "")), d))
  summary(m2 <- lm(as.formula(paste("ptv_SP.A ~ lotteryvariance.SP.A + lotterymean.SP.A + rat.SP.A + sex + as.factor(edu) + age", 
                                    if (robustnesscheck_pid) "+ pid_SP.A" else "")), d))
  summary(m3 <- lm(as.formula(paste("ptv_CDV ~ lotteryvariance.CDV + lotterymean.CDV + rat.CDV + sex + as.factor(edu) + age", 
                                    if (robustnesscheck_pid) "+ pid_CDV" else "")), d))
  summary(m4 <- lm(as.formula(paste("ptv_CDH ~ lotteryvariance.CDH + lotterymean.CDH + rat.CDH + sex + as.factor(edu) + age", 
                                    if (robustnesscheck_pid) "+ pid_CDH" else "")), d))
  summary(m5 <- lm(as.formula(paste("ptv_MR ~ lotteryvariance.MR + lotterymean.MR + rat.MR + sex + as.factor(edu) + age", 
                                    if (robustnesscheck_pid) "+ pid_MR" else "")), d))
  summary(m6 <- lm(as.formula(paste("ptv_OpenVLD ~ lotteryvariance.OpenVLD + lotterymean.OpenVLD + rat.OpenVLD + sex + as.factor(edu) + age", 
                                    if (robustnesscheck_pid) "+ pid_OpenVLD" else "")), d))
  summary(m7 <- lm(as.formula(paste("ptv_Ecolo ~ lotteryvariance.Ecolo + lotterymean.Ecolo + rat.Ecolo + sex + as.factor(edu) + age", 
                                    if (robustnesscheck_pid) "+ pid_Ecolo" else "")), d))
  summary(m8 <- lm(as.formula(paste("ptv_Groen ~ lotteryvariance.Groen + lotterymean.Groen + rat.Groen + sex + as.factor(edu) + age", 
                                    if (robustnesscheck_pid) "+ pid_Groen" else "")), d))
  
  BEL_2014_BNES_PS_Variance_Estimate <- tidy(m1) %>% filter(term=="lotteryvariance.PS") %>% select("estimate")
  BEL_2014_BNES_PS_Variance_SE <- se_robust(m1)["lotteryvariance.PS"] #tidy(m1) %>% filter(term=="lotteryvariance.PS") %>% select("std.error")
  BEL_2014_BNES_PS_Mean_Estimate <- tidy(m1) %>% filter(term=="lotterymean.PS") %>% select("estimate")
  BEL_2014_BNES_PS_Mean_SE <- se_robust(m1)["lotterymean.PS"] #tidy(m1) %>% filter(term=="lotterymean.PS") %>% select("std.error")
  
  BEL_2014_BNES_SP.A_Variance_Estimate <- tidy(m2) %>% filter(term=="lotteryvariance.SP.A") %>% select("estimate")
  BEL_2014_BNES_SP.A_Variance_SE <- se_robust(m2)["lotteryvariance.SP.A"] #tidy(m2) %>% filter(term=="lotteryvariance.SP.A") %>% select("std.error")
  BEL_2014_BNES_SP.A_Mean_Estimate <- tidy(m2) %>% filter(term=="lotterymean.SP.A") %>% select("estimate")
  BEL_2014_BNES_SP.A_Mean_SE <- se_robust(m2)["lotterymean.SP.A"] #tidy(m2) %>% filter(term=="lotterymean.SP.A") %>% select("std.error")
  
  BEL_2014_BNES_CDV_Variance_Estimate <- tidy(m3) %>% filter(term=="lotteryvariance.CDV") %>% select("estimate")
  BEL_2014_BNES_CDV_Variance_SE <- se_robust(m3)["lotteryvariance.CDV"] #tidy(m3) %>% filter(term=="lotteryvariance.CDV") %>% select("std.error")
  BEL_2014_BNES_CDV_Mean_Estimate <- tidy(m3) %>% filter(term=="lotterymean.CDV") %>% select("estimate")
  BEL_2014_BNES_CDV_Mean_SE <- se_robust(m3)["lotterymean.CDV"] #tidy(m3) %>% filter(term=="lotterymean.CDV") %>% select("std.error")
  
  BEL_2014_BNES_CDH_Variance_Estimate <- tidy(m4) %>% filter(term=="lotteryvariance.CDH") %>% select("estimate")
  BEL_2014_BNES_CDH_Variance_SE <- se_robust(m4)["lotteryvariance.CDH"] #tidy(m4) %>% filter(term=="lotteryvariance.CDH") %>% select("std.error")
  BEL_2014_BNES_CDH_Mean_Estimate <- tidy(m4) %>% filter(term=="lotterymean.CDH") %>% select("estimate")
  BEL_2014_BNES_CDH_Mean_SE <- se_robust(m4)["lotterymean.CDH"] #tidy(m4) %>% filter(term=="lotterymean.CDH") %>% select("std.error")
  
  BEL_2014_BNES_MR_Variance_Estimate <- tidy(m5) %>% filter(term=="lotteryvariance.MR") %>% select("estimate")
  BEL_2014_BNES_MR_Variance_SE <- se_robust(m5)["lotteryvariance.MR"] #tidy(m5) %>% filter(term=="lotteryvariance.MR") %>% select("std.error")
  BEL_2014_BNES_MR_Mean_Estimate <- tidy(m5) %>% filter(term=="lotterymean.MR") %>% select("estimate")
  BEL_2014_BNES_MR_Mean_SE <- se_robust(m5)["lotterymean.MR"] #tidy(m5) %>% filter(term=="lotterymean.MR") %>% select("std.error")
  
  BEL_2014_BNES_OpenVLD_Variance_Estimate <- tidy(m6) %>% filter(term=="lotteryvariance.OpenVLD") %>% select("estimate")
  BEL_2014_BNES_OpenVLD_Variance_SE <- se_robust(m6)["lotteryvariance.OpenVLD"] #tidy(m6) %>% filter(term=="lotteryvariance.OpenVLD") %>% select("std.error")
  BEL_2014_BNES_OpenVLD_Mean_Estimate <- tidy(m6) %>% filter(term=="lotterymean.OpenVLD") %>% select("estimate")
  BEL_2014_BNES_OpenVLD_Mean_SE <- se_robust(m6)["lotterymean.OpenVLD"] #tidy(m6) %>% filter(term=="lotterymean.OpenVLD") %>% select("std.error")
  
  BEL_2014_BNES_Ecolo_Variance_Estimate <- tidy(m7) %>% filter(term=="lotteryvariance.Ecolo") %>% select("estimate")
  BEL_2014_BNES_Ecolo_Variance_SE <- se_robust(m7)["lotteryvariance.Ecolo"] #tidy(m7) %>% filter(term=="lotteryvariance.Ecolo") %>% select("std.error")
  BEL_2014_BNES_Ecolo_Mean_Estimate <- tidy(m7) %>% filter(term=="lotterymean.Ecolo") %>% select("estimate")
  BEL_2014_BNES_Ecolo_Mean_SE <- se_robust(m7)["lotterymean.Ecolo"] #tidy(m7) %>% filter(term=="lotterymean.Ecolo") %>% select("std.error")
  
  BEL_2014_BNES_Groen_Variance_Estimate <- tidy(m8) %>% filter(term=="lotteryvariance.Groen") %>% select("estimate")
  BEL_2014_BNES_Groen_Variance_SE <- se_robust(m8)["lotteryvariance.Groen"] #tidy(m8) %>% filter(term=="lotteryvariance.Groen") %>% select("std.error")
  BEL_2014_BNES_Groen_Mean_Estimate <- tidy(m8) %>% filter(term=="lotterymean.Groen") %>% select("estimate")
  BEL_2014_BNES_Groen_Mean_SE <- se_robust(m8)["lotterymean.Groen"] #tidy(m8) %>% filter(term=="lotterymean.Groen") %>% select("std.error")
  
  # Harmonize names of the IVs of interest in the models
  names(m1$coefficients)[names(m1$coefficients) == "lotteryvariance.PS"] <- "Government Lottery Variance"
  names(m1$coefficients)[names(m1$coefficients) == "lotterymean.PS"] <- "Government Lottery Mean"
  names(m2$coefficients)[names(m2$coefficients) == "lotteryvariance.SP.A"] <- "Government Lottery Variance"
  names(m2$coefficients)[names(m2$coefficients) == "lotterymean.SP.A"] <- "Government Lottery Mean"
  names(m3$coefficients)[names(m3$coefficients) == "lotteryvariance.CDV"] <- "Government Lottery Variance"
  names(m3$coefficients)[names(m3$coefficients) == "lotterymean.CDV"] <- "Government Lottery Mean"
  names(m4$coefficients)[names(m4$coefficients) == "lotteryvariance.CDH"] <- "Government Lottery Variance"
  names(m4$coefficients)[names(m4$coefficients) == "lotterymean.CDH"] <- "Government Lottery Mean"
  names(m5$coefficients)[names(m5$coefficients) == "lotteryvariance.MR"] <- "Government Lottery Variance"
  names(m5$coefficients)[names(m5$coefficients) == "lotterymean.MR"] <- "Government Lottery Mean"
  names(m6$coefficients)[names(m6$coefficients) == "lotteryvariance.OpenVLD"] <- "Government Lottery Variance"
  names(m6$coefficients)[names(m6$coefficients) == "lotterymean.OpenVLD"] <- "Government Lottery Mean"
  names(m7$coefficients)[names(m7$coefficients) == "lotteryvariance.Ecolo"] <- "Government Lottery Variance"
  names(m7$coefficients)[names(m7$coefficients) == "lotterymean.Ecolo"] <- "Government Lottery Mean"
  names(m8$coefficients)[names(m8$coefficients) == "lotteryvariance.Groen"] <- "Government Lottery Variance"
  names(m8$coefficients)[names(m8$coefficients) == "lotterymean.Groen"] <- "Government Lottery Mean"
  
  
  if(!robustnesscheck_pid){
    
    stargazer(m1, m2, m3, m4, m5, m6, m7, m8, title="Linear regressions of vote choice on perceived government lottery variance and mean (Belgian National Election Study 2014)", 
              align=TRUE, 
              dep.var.labels=c("PS", "SP.A", "CD&V", "CDH", "MR", "OpenVLD", "Ecolo", "Groen"),
              omit.stat=c("LL","ser","f"),
              se = lapply(list(m1, m2, m3, m4, m5, m6, m7, m8), se_robust),
              p = lapply(list(m1, m2, m3, m4, m5, m6, m7, m8), p_robust),
              out = "TableSM6.tex"
    )
    
  }

